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INVERSE SOLUTIONS FOR LAMINAR BOUNDARY- LAYER FLOWS 


WITH SEPARATION AND REATTACHMENT 

James E. Carter 
Langley Research Center 

SUMMARY 

Numerical solutions of the laminar, incompressible boundary -layer equations are 
presented for flows involving separation and reattachment. Regular solutions are obtained 
with an inverse approach in which either the displacement thickness or the skin friction is 
specified; the pressure is deduced from the solution. A vorticity — stream -function formu- 
lation of the boundary -layer equations is used to eliminate the unknown pressure. Solu- 
tions of the resulting finite-difference equations, in which the flow direction is taken into 
account, are obtained by several global iteration schemes which are stable and have 
unconditional diagonal dominance. Results are compared with KJineberg and Steger's 
separated boundary-layer calculations, and with Briley's solution of the Navier-Stokes 
equations for a separated region. In addition, an approximate technique is presented in 
which the streamwise convection of vorticity is set equal to zero in the reversed flow 
region; such a technique results in a quick forward- marching procedure for separated 
flows. 


INTRODUCTION 

Despite the recent advances made in computer technology and computational tech- 
niques, numerical solutions of the Navier-Stokes equations for separated flows are still 
impractical for many applications and in many instances are unnecessary. There is 
growing evidence that the boundary-layer equations provide a reasonably accurate model 
for separated flows of limited extent. For example, good agreement with experiment was 
obtained by Klineberg and Lees (ref. 1) by using an integral technique to solve the 
boundary-layer equations for a supersonic, separated flow over a compression corner. 
More recently, similar calculations were made by Dwoyer (ref. 2) and by Werle and 
Vatsa (ref. 3) using finite-difference techniques which were in good agreenient with exper- 
imental results as well as with the Navier-Stokes computations made by Carter (ref. 4). 
Similarly, Ghia and Davis (ref. 5) concluded from their numerical solutions of the incom- 
pressible Navier-Stokes equations for separated flows that the boundary-layer equations 
appear to be adequate for flows with small separation regions provided that displacement 



effects are appropriately considered. These results warrant the further examination of 
the applicability of the boundary-layer equations for describing separated flows. In the 
present report some additional understanding of their applicability is provided through 
the examination of results obtained by two solution techniques for the laminar boundary- 
layer equations for incompressible flow. 

As pointed out by Brown and Stewartson (ref. 6), it is generally agreed that the 
solution of the boundary- layer equations with a prescribed pressure gradient results in a 
singularity at the separation point. Many numerical investigations such as those made by 
Terrill (ref. 7), by Werle and Davis (ref. 8), and by Klineberg and Steger (ref. 9) support 
this conclusion. In all of these investigations, good agreement was obtained with 
Goldstein's (ref. 10) analytical results near separation. In supersonic flow the singular- 
ity at separation can be eliminated by modifying the pressure distribution through viscous - 
inviscid interaction. At each streamwise station the pressure calculation is coupled to 
the growth of the boundary layer through the use of the Prandtl- Meyer or tangent-wedge 
relationship. In subsonic flow the inviscid flow is elliptic. Therefore, such a procedure 
is not possible because the surface-pressure distribution is dependent on the entire body 
shape; the entire boundary-layer flow has to be taken into account before its effect on the 
inviscid flow can be calculated. 

Several numerical investigations have demonstrated that regular solutions of the 
boundary -layer equations at separation can be obtained (excluding interaction) by specify- 
ing either the displacement thickness or the wall-shear stress distribution; the pressure 
distribution is deduced from the resulting solution. Catherall and Mangier (ref. 11) were 
the first to demonstrate that the boundary -layer equations could be numerically integrated 
past the separation point (without a singularity) by prescribing the displacement thickness 
distribution. With the wall shear prescribed, Kuhn and Nielsen (ref. 12) presented calcu- 
lations in which they used an integral technique to solve the turbulent boundary -layer 
equations for separated flow. Klineberg and Steger (ref. 9) and, more recently, Horton 
(ref. 13) used finite-difference techniques to solve the laminar boundary -layer equations 
with the wall shear prescribed. In all of these calculations the singularity usually pres- 
ent at separation was eliminated. Another inverse procedure, which was developed ear- 
lier by Klineberg and Steger (ref. 14) and later used by Tai (ref. 15), is to use the trans- 
verse velocity at the boundary- layer edge as the prescribed condition. Since the edge 
value of the normal velocity component is related to the streamwise growth of the dis- 
placement thickness, this procedure is somewhat analogous to specifying the displacement 
thickness. 

Inverse procedures have been used for attached flow as well as for separated flow. 
For example, Keller and Cebeci (ref. 16) developed an inverse calculation procedure for 
calculating attached laminar flow for a specified wall shear. Cebeci and Witherspoon 
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(ref. 17) modified the procedure to study incipient turbulent separation including the 
effects of suction at the wall. 

In contrast with these inverse techniques, Briley and McDonald (ref. 18) have made 
calculations for subsonic flow using a glirect procedure. In this procedure the unsteady 
boundary -layer equations are repeatedly solved until a steady-state solution is obtained. 
After each time step, the prescribed pressure is updated from thin airfoil theory, thereby 
accounting for the displacement thickness interaction. Although this technique seems 
feasible, it needs further examination since Briley and McDonald obtained a regular solu- 
tion at a laminar separation point for a case with no interaction. This result is probably 
because of their use of a first-order scheme in the streamwise convection terms. As 
shown by Werle and Davis (ref. 8) and Klineberg and Steger (ref. 9), such a scheme 
requires a very fine mesh in order to resolve the separation singularity. Hence, in those 
cases in which interaction was included, it is not clear whether the solution at separation 
would be regular if a second-order scheme were used. 

In the present report, which is based in part on references 19 and 20, numerical 
solutions of the incompressible, laminar boundary-layer equations for flows with separa- 
tion and reattachment are presented. Solution techniques for either a prescribed dis- 
placement thickness or wall shear are described. In a complete calculation of a viscous- 
inviscid interaction, this boundary -layer computation would be matched to a calculation of 
the inviscid flow and this process then iterated until convergence. In the present report 
no viscous-inviscid interaction is considered; however, in a recent paper by Carter and 
Wornom (ref. 21), the prescribed displacement thickness method was combined with an 
inverse inviscid-flow analysis to compute a subsonic viscous-inviscid flow interaction. 

In addition to the interaction considerations, prescribing the displacement thickness pro- 
vides a useful scaling function for the normal coordinate y, thereby eliminating boundary- 
layer growth in the computational coordinate r] = y/5*, where 6* is the displacement 
thickness and r] is the transformed normal coordinate. Therefore, this approach has 
received the main emphasis in the present study, although some calculations, mainly for 
comparative purposes, have been made with the wall shear prescribed. With the wall 
shear prescribed, several calculations were found to give discontinuous solutions. In an 
effort to explain these discontinuities, three different formulations of the boundary -layer 
equations for a prescribed wall shear are considered. 

Gather all and Mangier (ref. 11), in their numerical study of separated flow, did not 
account for the reversed flow direction in their finite-difference scheme. In reference 11, 
the same finite-difference scheme was used in the reversed flow region as the scheme 
used for attached flow; hence, only one streamwise pass was made through the boundary 
layer. Since integration of the boundary -layer equations may be unstable in a direction 
opposite to that of the flow, it is not surprising that Catherall and Mangier encountered 
problems in convergence of the column iteration in the reversed flow region. In the 
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present report the lack of convergence found by Catherall and Mangier in the reversed flow 
region is overcome either by the use of a global iterative procedure or by an approximate 
forward-marching technique. In the global scheme, the finite-difference approximation for 
the streamwise convection term is switched from a backward to a forward difference in the 
reversed flow region to account properly for the flow direction. Klineberg and Steger 
(ref. 9) used a similar procedure in their numerical computation of separated flow. With 
the global scheme, the solution region must extend to a point downstream of reattachment 
to avoid the requirement of specifying the unknown downstream boundary conditions in the 
reverse flow region. This scheme requires repeated sweeps from the upstream boundary 
through reattachment until convergence is obtained. 

An approximate forward- marching scheme is presented for a prescribed displace- 
ment thickness in which the streamwise convection of vorticity is neglected in the reverse 
flow region. The resulting procedure for obtaining approximate solutions for separated 
flows is quite rapid and requires only modest computer storage since the procedure is 
essentially the same as the usual forward- marching technique used to solve attached 
flows. This approximation is similar to that used by Reyhner and Fliigge-Lotz (ref. 22); 
they were the first to show that a stable calculation can be obtained in marching through a 
reversed flow region if the streamwise convection of momentum is set equal to zero. 

Werle, Polak, and Bertke (ref. 23) and Williams (ref. 24) have also used this approxima- 
tion in their calculations. 

This report is divided into a number of major headings and subsections to aid the 
reader. In the first two major headings after the symbols list, the prescribed displace- 
ment thickness and the prescribed wall-shear methods are discussed, respectively. In 
these sections both the mathematical formulations as well as the numerical procedures 
are presented. The stability of the various finite -difference schemes used in these 
inverse methods is examined with the von Neumann analysis; this analysis is presented 
in appendixes A and B. The results obtained with these two inverse procedures are then 
presented separately in the "Results and Discussion" section. 

SYMBOLS 


An»Bnri 
Cn>Dn J 


coefficients in tridiagonal system of equations 


Cf skin-friction coefficient 

Cn,Dn coefficients in Thomas algorithm 

Cz coefficient defined in equation (81) 
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coefficient defined in equation (20) 




f 


g 

J 

L 


M 


CO 


m 


m ,n 


q 

^°o,L 


r 

S 

t 

At 

Uoc 

u 

u 

% 


coefficient defined in equation (19) 

coefficients in recurrence equation for stream function 

scaling function for normal coordinate 

amplification factor 

global iteration counter 

reference length 

free-stream Mach number 

pressure-gradient parameter 

indices for and r]-directions, respectively 

column iteration counter 

free-stream Reynolds number, 

^ CO 

relaxation factor 

stream function obtained from column solution 
time 

time increment 
free-stream velocity 

velocity component parallel to surface divided by Uoo 

_ _u_ 

" Ue 

boundary velocity divided by Uoo 


velocity component parallel to surface at boundary -layer edge divided by Uoo 


velocity component normal to surface divided by Uoo 
= 

transformed normal velocity component 
vorticity obtained from tridiagonal equations 
coordinate along surface divided by L 
coordinate normal to surface divided by L 

= ypoo^L 

transformed normal coordinate 
coefficient of artificial time term 
parameter in prescribed wall shear 
displacement thickness divided by L 



transformed normal coordinate 
grid spacing in 77-direction 
^ argument in von Neumann analysis 
Meksyn pressure -gradient parameter 
molecular viscosity coefficient 
transformed x- coordinate 


grid spacing in | -direction 



free- stream density 


Poo 

f transformed shear 

Tq transformed wall shear 

(f) 7] argument in von Neumann analysis 

\l/ stream fvmction divided by UooL 

transformed stream function in prescribed displacement thickness method 
}p transformed stream function in prescribed wall- shear method 

(jO vorticity divided by Uco/L 

cJq prescribed wall vorticity divided by Uoo/L 

PRESCRIBED DISPLACEMENT THICKNESS METHOD 

Governing Equations and Boundary Conditions 

The nondimensional boundary-layer equations for two-dimensional, laminar, incom- 
pressible flow can be written as 


Bu ^ 9y 
9x ay 


= 0 



- dxi 

+ V = Up 

9y 




with the boundary conditions 


u(x,0) = v(x,0) = 0 

and 

u(x,y) - Ug(x) as y - 


( 1 ) 

( 2 ) 


(3) 

(4) 
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where Ug is the velocity at the edge of the boundary layer. All barred quantities have 
been scaled in the usual manner by the Reynolds number and are given in detail in the 
list of symbols. The detailed derivation and discussion of these equations and boundary 
conditions are presented by SchUchting (ref. 25). 

This formulation is modified for the present application by first eliminating the 
unknown edge velocity Ug from the problem. Taking the y-derivative equation (2) and 
introducing the vorticity To = 9u/9y results in 


u + V 

9y 8y2 


(5) 


In order to eliminate the edge velocity from the outer boundary condition, it is convenient 
to introduce the stream function 





( 6 ) 



J 

The first part of equation (6) can be combined with the displacement thickness 



to show that as y 

^ - Ug(y - 5*) (8) 

Therefore, if a transformed stream function is introduced 


= ^ - u(y - 6*) (9) 

then equation (4) may be replaced by 

1 ^ — 0 as y — (10) 
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If the independent variables are transformed by 





and if equations (6) and (9) are combined and substituted into equations (1) and (2), then the 
governing equations become 



9a; 



+ (?7 - l)(u6*) 


8u) 

dj] 



(12) 


=5* (1 - T])u) 
dri 

subject to the boundary conditions 


(13) 


uU,0) = ^{^,0) = 0 


(14) 


w(|,77) - 0 

M^,v) - 0 


as T 7 — °o 


(15) 


The transformed stream function given in equation (9) does not eliminate the u compo- 
nent of velocity from the governing equations. However, once the vorticity is deduced 
from equation (12), then u is given by 


u(^,T?) = 6* w(^,?7l) d?7i (16) 

where u(^,0) = 0 has been used. For a given displacement thickness distribution, the 
use of equations (12), (13), and (16) subject to the boundary conditions of equations (14) 
and (15) results in a description of the boundary-layer flow in which the explicit appear- 
ance of the unknown edge velocity Ug has been eliminated. Unless otherwise noted, in 
the remainder of this paper the bars, which have been previously used to denote Reynolds 
number scaling, are deleted. 
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Numerical Technique: Global Iteration 


Finite-difference scheme .- The governing equations are replaced by finite-difference 
expressions which depend on whether the flow is forward (u = 0) or reversed (u < 0). The 
presence of the reverse flow region requires repeated sweeps from the upstream bound- 
ary, which is located upstream of the separation point, to a point downstream of reattach- 
ment. This sweep procedure is termed a global iteration scheme and is described in this 
section. Figure 1 gives a schematic diagram which shows the computational molecules 
used in the forward and reversed flow regions for the vorticity equation both for the pre- 
scribed displacement thickness approach and for one of several wall -shear prescribed 
approaches. The figure also specifies the boundary conditions. The indices m and n 
denote ^ = m and 77 = n A?7 and the cross in the center of the computational mole- 
cule denotes the point at which equation (12) is evaluated. In the forward flow regions it 
was found to be imperative to maintain second-order accuracy, and thus equation (12) is 

evaluated with the Crank-Nicholson scheme at the point | = ^m - 77 = n A77. 


Forward flow Reversed flow 


Displacement 

thickness 

prescribed 

(DTP) 




Prescribed 

upstream 

flow 


Wall shear 
prescribed 
(WSP) 




Outer boundary: co= 0, 0 (DTP only) 



No condition 
on downstream 
boundary 


Figure 1. - Computational schemes and boundary conditions used 
for boundary- layer equations. 
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In the reversed flow regions first-order accuracy in the ^ derivative was found to 
be adequate; hence, a first-order forward-difference expression is substituted for 
Instabilities were encountered when the same forward-difference approximation was used 
on and on 8u6*/8^, These instabilities are eliminated when central-difference 

expressions are used for these terms. In the reversed flow region, 8 oj/8t] and d^u>/Br]^ 
are approximated by central-difference expressions; the same approximation is used in the 
forward flow region. 

In both the forward and reversed flow regions the finite-difference representation of 
equation (12) can be written as 


■^n‘^m,n-l + ®n^m,n ^n^m,n+l " ^n 


If the flow is forward 

A„=-i(c,+ l) 

Bn ~ 1 + 2C ^ 

C„ = i(c, - 1) 

where 




. ^fu6*2' 
2 A^\ 



(17) 


(18) 


(19) 


C 


V ~ 



B 

H 


i// + (t 7 - l)(u6 




( 20 ) 


In equations (19) and (20) average values are used for the required quantities at 
I = ^m - In equation (20) a central difference is used for the | derivative. If the 

flow is reversed, then 
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An = -{prj + l) 

Bn = 2^1 - 

- ^T] ~ ^ 

Dn = ~2C^a)m+l^n 


where 





( 21 ) 


( 22 ) 


C 


V ~ 


a 

2 


^ + (v - 


l)(u5*) 


m,n 


(23) 


No blending between the forward and reversed flow difference expressions such as that 
used by Klineberg and Steger (ref. 9) was used in the present iterative calculations. 

Repeated application of equation (17) from the wall to the outer boundary results in 
a tridiagonal system of linear equations for the vorticity if the coefficients An, Bn, Cn, 
and Dn are assumed known. These equations are easily solved by the Thomas algorithm 
(ref. 26), which can be written as 

^m,n = d; + Cncojn n-i (24) 


where 


> (25) 

+ ^n^n+1 

y 

The quantities and are computed from the outer boundary to the wall with the 

outer boundary condition = 0 imposed by setting = 0, where N denotes 

the grid point at the outer boundary. Equation (24) is then used to deduce provided 

the value at the wall is known. An implicit procedure for deducing the wall vor- 

ticity is described in the next section. 


Dn - CnD: 


^n^n+1 
Dn + ^n^n+1 
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Calculation of wall vorticity .- The wall vorticity can either be estimated and later 
updated from the subsequent solution, or it can be found along with the column solution by 
simultaneously solving the stream-function equation (eq. (13)). The latter procedure was 
used in the present investigation. 

Equation (13) can be evaluated to second-order accuracy at the point ^ = m 

V = U + ^^V by 


'^m,n+l ” ^m,n 


Ar] 




^m,n+l + ^^m,n 
2 


(26) 


Equation (13) must be evaluated at a half-grid point above, and not below the level at which 
the vorticity equation (eq. (12)) is differenced so that the outer boundary condition 

= 0 enters the problem. Equation (26) can be combined with equation (24) to 
obtain a recurrence relation for the stream function which is given by 


’^m,n “ ^n + Ln^m,n-1 


(27) 


where 


^n “ ^^n+1 ^n 


d;+i + d;(i + c;+i) 


+ ^n^n+1 


Ln - Ln+1 + E^^l + 




- 1 


^^2 


J 


(28) 


The quantities and are computed from the outer boundary, to the wall, with the 

boundary conditions w(?,°°) = = 0 imposed by setting K]^ = Lj,j = 0. The wall 

vorticity j is found by combining equation (26), which is evaluated at rj = Arj/2 and 

in which the boundary condition i//(^,0) = 0 is imposed, with equations (24) and (27) which 
are evaluated at r\ - Arj (n = 2) to obtain 


f - K2 

1 ^ (29) 

^ + l) - L2 
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Using this result, equations (24) and (27) may now be used to compute ’^mjU 

by back substituting from the wall to the outer boundary. 

Stability considerations .- Several global iterative procedures for solving equa- 
tions (17) and (26) were developed, but were found to be unstable in certain calculations, 
particularly as the magnitude of the stream wise gradients increased. However, these 
instabilities were eliminated by the introduction of a time-like term in equation (17) to 
establish unconditional diagonal dominance. In addition, the introduction of this term 
results in a global iterative scheme which can be analyzed for linear stability by the von 
Neumann analysis. Furthermore, it was necessary to use underrelaxation to obtain con- 
verged solutions. Denoting the unknown vorticity in equation (17) as Wn, these modifi- 
cations can be described by rewriting equation (17) as 


^n^n-l (®n + “)^n + C^Wn+l - Dn + 
followed by the relaxation equation 


O)' 


J+1. = J 


m,n 


r W„ - wJ 


(30) 


(31) 


where the superscript j denotes the global iteration level and r is the relaxation fac- 
tor. Thus Wn is a temporary value of the vorticity along the column currently being 
updated. Equation (27) for the stream function is modified to give 


Sn = Kn + LnWn-i (32) 

with the relaxed value of the stream function found from 

= p + r(sn - P^^^ (33) 

Expressions for a for forward and reversed flow are determined so that diagonal 
dominance is maintained; that is. 


Bn + Q! 



"n 


(34) 


even if 


= 0, and these expressions are given by 


u s 0 
u < 0 





0 = 2 


C?7 






(35) 
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where is determined by equation (20) or (23) depending on whether the flow is for- 


ward or reversed. These values of a insure that if equation (34) is satisfied, then, as 
shown inductively by Keller (ref. 27), jc^l = 1, thereby insuring no growth of round-off 

errors in applying equation (24). In the present calculations, modified expressions for a 
were used such that | Ln | = 1 as well as j | = 1. These relationships are found to be 


u g 0 


u < 0 


a = C 




+ 1 1 + En 


ot — 2j|C7^ + 1 + C'q Ej^j 


(36) 


where, as before, is determined by either equation (20) or (23); is given in 

equation (28). 

The linear stability of the finite-difference forms of the vorti city -transport equa- 
tion is examined with the von Neumann analysis in appendix A. Both forward and reversed 
flow are considered. A stability analysis of the stream -function computation is considered 
unnecessary. The stream -function calculation, despite the complications introduced by the 
implicit treatment of the wall vorticity, is simply an integral of the vorticity across the 
boundary layer and should be stable (provided the vorticity computation is stable). The 
stability analysis of the vorticity computation is begun by eliminating the temporary value 
Wn between equations (30) and (31). The elimination results in 




J+1 






(37) 


In appendix A the details of the von Neumann analysis of equation (37) are presented. 

This analysis shows that as A|, A ?7 - 0, with At^2/a^ held constant to preserve the 
highest 4 3,nd rj derivatives, the only restriction for either forward or reversed flow 
is that r = 1. In this limit the lower order derivatives are suppressed by the zero grid 
size; the remaining terms form the relaxation solution of the diffusion equation 

(38) 

where a is the di^usion coefficient. Thus, with the present schemes the diffusion equa- 
tion cannot be overrelaxed. 
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For finite grid size the inclusion of the lower order terms has no effect on the first- 
order scheme used in reversed flow. However, for forward flow it is found that the lower 
order terms may be destabilizing, and thus may impose a limit on the maximum size of 
the grid spacings which can be used. This point is further discussed in appendix A along 
with a discussion of some numerical results which qualitatively confirm the stability 
analysis. However, it appears from the calculations that the grid size used in the pres- 
ent study is sufficiently small to suppress the adverse stability effects of the lower order 
terms. Hence, other than the restriction of under relaxation, the calculations presented in 
this report were performed with no grid constraints other than that of truncation errors. 
The effect of the lower order terms in the present investigation is the same as that dis- 
cussed by Richtmyer and Morton (ref. 28) in that, for diffusion problems, the stability 
criterion is essentially unaltered by the inclusion of the lower order terms for a wide 
range of finite-difference schemes. 

An alternate scheme for forward flow with second-order accuracy is also presented 
in appendix A. In this scheme, the destabilizing influence of the lower order terms is 
eliminated when no underrelaxation is used; that is, when r = 1. This scheme was de- 
scribed previously in reference 19 as having unconditional stability by the von Neumann 
analysis. This condition is now shown to be true only when no relaxation is used. How- 
ever, underrelaxation was generally used in the present calculations; thus, this alternate 
scheme is of little value. The value of the relaxation factor typically used in the present 
study is between 0.4 and 1.0. A short discussion based on some numerical results is 
given in appendix A in which it is suggested that underrelaxation is required primarily so 
that the computation of the surface vorticity is stabilized. 

Computation procedure .- The global iterative procedure is initialized by guessing 
the values of w and \p at each point in the computational region (fig. 1). At the up- 
stream boundary, the boundary conditions for w and xj/ are prescribed. At the down- 
stream boundary, it is assumed that the flow is attached; hence, it is unnecessary to 
impose any conditions here. The new iteration values of cu and ^ are computed from 
equations (31) and (33), respectively, along the column of grid points from the wall to the 
outer boundary. In this computation, the most recently updated values of cu, !^, and u 
have been used in evaluating the coefficients. The columns are computed successively 
beginning at the upstream boundary and continuing to some point downstream of reattach- 
ment. Since the flow is reversed and since an artificial time term has been introduced 
in the vorticity equation in order to maintain diagonal dominance, this sweep procedure 
from the upstream to the downstream boundary must be repeated iteratively. This 
global iterative process is continued until the following convergence criterion is satisfied: 


max 

m,n 



- f^ 
m,n 


^ 10"5 


(39) 
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where f represents co, rp, or u, that is, whichever quantity results in the largest 
residual. 


Forward- Marching Procedure 

In cases in which the velocity in the reversed flow region is small (5 percent or less 
of the free- stream velocity), the global iteration procedure described in the previous sec- 
tion is unnecessary. In these cases neglect of the streamwise convection in the reversed 
flow region should have only a slight influence on the resulting solution. For example, 
Reyhner and Fliigge-Lotz (ref. 22) demonstrated that, by setting the convection term 

u in the x momentum equation equal to zero for u less than zero, a stable finite- 

difference solution of the boundary-layer equations could be obtained with the usual 
forward- marching procedure for a separated flow. This approximation eliminates the 
well-known instability encountered in solving the boundary-layer equations in a direction 
opposite to that of the local flow. The approximation also results in a substantial reduc- 
tion in computer time and storage as compared to that required for a global iteration pro- 
cedure. Several investigators (Dwoyer (ref. 2), Werle, Polak, and Berthe (ref. 23), and, 
more recently, Williams (ref. 24)) have used the Reyhner and Fliigge-Lotz approximation 
in their separated boundary -layer calculations. An approximation similar to that of 
Reyhner and Fliigge-Lotz (discussed by Carter and Wornom in ref. 20) is made in the 
present formulation by neglecting the streamwise convection of vorticity in the reversed 

flow region. Thus, when u < 0, set u6*^ — = 0; this approximation is implemented by 
setting C| = 0. 

With the inclusion of this approximation, the global iteration technique becomes the 
usual forward- marching procedure with column iteration which is typically used to solve 
the boundary- layer equations for an attached flow. Instabilities encountered in these 
calculations were eliminated by introducing time-like terms in the finite-difference equa- 
tions, similar to those terms used in the global procedure. The introduction of these 
terms provides both unconditional diagonal dominance and a column iterative scheme which 
was found to be stable using the von Neumann stability analysis. In some cases Reyhner 
and Fliigge-Lotz (also Werle, Polak, and Bertke (ref. 23) in a similar investigation) 
encountered an unexplained instability in the reverse flow region; they eliminated this 
instability by introducing a positive artificial convection term. Since this term increases 
the magnitude of the diagonal coefficients in the tridiagonal system of equations, it is pos- 
sible that the instability in their solutions was caused by large errors in the Thomas algo- 
rithm because of a lack of diagonal dominance. 
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The approximate forward- marching technique with column iteration is somewhat 
similar to the global iterative procedure. Unconditional diagonal dominance is estab- 
lished by introducing a time-like term in equation (17) to give 

An^^Wn-i + (Bn + a)%n + Cn%n+1 = + (a‘^m,n)^ (40) 


where q denotes the column iteration level. The value of a which was used in the 
calculations is given by equation (36) for u = 0. Underrelaxation was used between re- 
peated column iterations and is written as 


m,n m,n 


+ r(Wn 


(jO 


q 

'm,n 


(41) 


This modification of the finite-difference equations to obtain diagonal dominance 
results in a column iterative procedure which may be approximately analyzed by the von 
Neumann analysis. As in the global iterative procedure, the temporary value is 

eliminated between equations (40) and (41) to give 




m,n- 


(B„ H- = rDn'l . . (1 - r)( ^ 




+ B 


^ 0)^1 + C „ 

m,n n 


n 


‘leu? 


m,n 


(42) 


The substitution of a single Fourier component of arbitrary wave number k, 


cu 




n 


into equation (42), where cuq is the amplitude, results in the amplification factor 


cu^^l 

0 


Q! + (1 

- r)^l - cos 0 + + i(l - r)Cj^ sin 4> 



0! + 1 - COS 0 + 2C^ + iCTj sin 0 



(43) 


where <p = yi Arj and i = \pl. The derivation of equation (43) is simplified by assuming 
Dn‘1 is a constant. Since Dn*l depends on A^, Cn, and so on, which are assumed to be 
constant, it is neglected in the stability analysis. The von Neumann condition for stability 
is that jgj § 1. Since =0 and a = 0, this condition is satisfied in equation (43) if 
r = 1. Thus, the repeated column iterations have linear stability provided that over relax- 
ation is not used; the same result was obtained previously for the global iterative proce- 
dure. This technique for obtaining unconditional diagonal dominance and column stability 
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is general, and should be applicable to other implicit finite-difference schemes such as 
the alternating direction implicit (ADI) solution procedure for solving the Navier-Stokes 
equations. 


PRESCRIBED WALL-SHEAR FORMULATION 


Governing Equations and Boundary Conditions 

Three formulations of the boundary -layer equations for a prescribed wall shear are 
given below. The first formulation was given in reference 19; in certain calculations it 
resulted in a solution which varied discontinuously in the x-direction at a particular x 
location. More recently it has been found that this discontinuity can be eliminated by 
introducing a constraint on the edge velocity Ug, which is described in the second 
formulation. 

The third formulation is a vorticity — stream-function equivalent of Klineberg and 
Steger's (ref. 9) formulation of the boundary-layer equations. The third formulation is 
expressed in primitive variables u and v. In one of Klineberg and Steger's calcula- 
tions, an unexplained discontinuity occurred in the pressure-gradient parameter 
^ X ^^e 

m = which appears explicitly in their governing equations. In reference 19 it was 

conjectured that the two discontinuities mentioned might be related. This report shows 
that there is no relationship; hence, the discontinuity encountered by Klineberg and Steger 
remains unexplained. In an attempt to understand this discontinuity further, this case was 
recomputed using the present solution technique. However, in order to specify the same 
transformed wall-shear distribution as Klineberg and Steger, it was necessary to develop 
a new vorticity — stream-function formulation (designated as formulation 3). 

Formulation 1 (co, u) .- The boundary -layer problem with the wall shear pre- 

scribed can be formulated in an analogous manner to that with the displacement thickness 
prescribed. In contrast to the primitive variables used by Klineberg and Steger (ref. 9), 
the use of vorticity as a dependent variable seems particularly simple for a prescribed 
wall shear since Dirichlet boundary conditions on the vorticity are imposed at both the 
wall and boundary-layer edge. The governing equations and boundary conditions for a 
prescribed wall- shear distribution Wq(^) can be written as follows, where the displace- 
ment thickness (now unknown) has been replaced in equations (12) and (13) by a given 
scaling function f(x), 


uf 


2 9a) 




\j/ + (t] - 1) (uf) 


8cu 


dri^ 


(44) 
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(45) 


B\l/ 

dj] 


= f^(l - ri)oi) 


u(|,0) = ^{^,0) = 0 
o)(4,0) = o)o(^) 


(46) 


aj(|,77)-0 as 

where the independent variables are 



and the transformed stream function is 


(47) 


(48) 


4/ = \l/ - uf(7] - 1) 


(49) 


After the vorticity is deduced, u is given by 


u(|,7?) = f yj <^U,Vi) dVi (50) 

where the boundary condition u(^,0) = 0 has been imposed. In the present calculations 
the scaling function is chosen as f(x) = 6 q*^x/xq where 6 q* is the displacement thick- 
ness at the upstream boundary located at Xq. 

After the calculation is completed, the unknown edge velocity may be deduced from 
the velocity profile given by equation (50) with 77 = °° or from the x momentum equa- 
tion evaluated at 77 = 0 which gives 


u 


e 


dUe 


1 8gj 

f 977 


'77=0 


(51) 


The edge velocity is deduced from equation (51) by integrating from the upstream bound- 
ary where Ug is known. The value of Ug (referred to as the profile value) deduced 
from equation (50) and that value deduced from equation (51) should agree; in fact, this 
redundancy serves as a useful check on the solution accuracy. For those solutions in 
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which discontinuities occurred, it was found that the Ue from the profile was discon- 
tinuous; that value deduced from equation (51) was not. Therefore, this discontinuity was 
eliminated by forcing these two Ue values to be equal. This step results in the modified 
formulation presented in the following section. 

Formulation 2 (co, u; Ue).- The transformed stream function is replaced with 
the usual stream function in equation (44), and the vorticity-transport equation is then 
written as 


du) dip 8cu\ _ 

\8t7 

subject to the boundary conditions 

w(|,0) = cooiO 

~ 0 as T 7 — 


(52) 


(53) 

(54) 


After the vorticity is deduced, the u component of velocity is found from the solution of 
= (55) 

Bri^ 

Equation (55) is obtained by differentiating the definition of vorticity w =i— with re- 

I Brj 

spect to 77 . It is convenient to solve this second-order equation for the velocity in order 
to impose simultaneously both of the following boundary conditions: 


and 


u(?,0) = 0 


u(^,?7) - ue(l) = 



1I/2 




as rj 


(56) 


(57) 


The latter of these conditions is the integral of equation (51) where is the position of 
the upstream boundary. The stream function is then given by 


rv 

= f ii(^,7?l) 

where the boundary condition i/^(|,0) = 0 has been imposed. 


(58) 
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In order for equation (57) to be satisfied, the vorticity must vanish at large distances 
from the surface. Thus, the imposition of equation (57) automatically insures that equa- 
tion (54) is satisfied. Hence, the main difference between the two formulations given 
herein is that, in the second, a stronger outer boundary condition is imposed by constrain- 
ing Ue as well as imposing a zero value of vorticity at the boundary -layer edge. In the 
first formulation, only the latter boundary condition is imposed. 

Formulation 3 (t, ij/, u; m ).- Klineberg and Steger (ref. 9) solved the following 
boundary -layer equations 


0X dy 2 


and 


where 




= m(l - u2) 




(59) 


(60) 


(61) 


Klineberg and Steger solved these equations subject to the boundary conditions 
y = 0 u = v = 0 = Tg(x) 

9y 


y 


CO 


u - 1 


(62) 

(63) 


where Tq(x) is a prescribed function. The pressure-gradient parameter m is deduced 
from equation (60) evaluated at y = 0; this step gives 


m = 



(64) 
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The transformed shear r is related to usual vorticity co by 


T = 



( 65 ) 


which involves the unknown edge velocity. Therefore, in order to compute the same case 
as Klineberg and Steger (ref. 9) with the present technique, it is necessary to write the 
vorticity equation in terms of the transformed shear r. This step is easily done by dif- 
ferentiating equation (60) with respect to y. The differentiation gives 


d\j/ dr 

3y 9x 




.2- 


8t ^ _ - 8ii/ /3m - 1 

9y 9 A 2 / 


( 66 ) 


where the stream function has been introduced as 


^ Bii/ 

u = — ^ 

9y 


v = -Al + m)-x|^ 
2 9x 




(67) 


Equation (67) automatically insures that the continuity equation is satisfied. The trans- 
formed shear is determined from equation (66) subject to the boundary conditions 


t(x,0) = Tq(x) 


(68) 


T(x,y) - 0 as y - «> 

(69) 

The remainder of the formulation closely follows that given in formulation 2. 
transformed shear is deduced, the u component of velocity is found from 

After the 

9^u _ ^ 

C70) 

3y2 9y 


subject to the boundary conditions 


u(x,0) = 0 

(71) 

u(x,y) - 1 as y - " 

(72) 
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The stream function is then given by 


$^(x,y) = u(x,y J dyj 


( 73 ) 


where the boundary condition i^(x,0) = 0 has been imposed. The pressure-gradient 
parameter is then found from 


m(x) = 

9y 


y=0 


( 74 ) 


Numerical Technique 

Formulation 1 .- The finite-difference form of equations (44) and (45) is solved with 
a global iteration procedure similar to that presented for the prescribed displacement 
thickness formulation. However, in the forward flow regions, the Crank-Nicholson 
scheme produced small streamwise oscillations which were eliminated by use of the com- 
putational molecule shown in figure 1. In this case the ^ derivatives are represented 
to second-order accuracy by a three-point backward-difference expression; the rj deriv- 
atives in the vorticity equation are represented by the usual central-difference expres- 
sions. In the reversed flow region, the finite-difference representation of the vorticity 
equation is the same as that used for the prescribed displacement thickness approach. 

In both the forward and reversed flow regions, the finite-difference representation 
of equation (44) can be written in the form of equation (30) where, if the flow is forward 


■^n - ~[^ri + 1 


+ Of — 2 + 3Ct + 2 


^n “ ^77 ^ 


( 75 ) 


and 


Dn + QfwJ „ = Ct(4a)J+\ „ - „ ) + 2 

“ m,n m-ljii m-Z,nJ 


Ct)-’ 


C,=J^uf2 


• = ^lUI / 

^ 2A^^ ^m,n 


~2Xh' 


rp + in - l)(uf) 




( 76 ) 
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Three-point backward-difference expressions are used for the | derivatives in equa- 
tion (76). If the flow is reversed 




An = 

-(c» ^ 1) 


Bn + 

Q! = 2^1 - 



o 

+ 

u 

Cn = 

C„- 1 


Dn + 

“^m,n = 





and 

C-q are given in equation (76). For 


(77) 


are used for the 4 derivatives in C 


V 


The solution of equation (30) for the vorticity is again found with the Thomas algor- 
ithm. The Thomas algorithm is given in equations (24) and (25) subject to the boundary 
conditions of equations (46) and (47). The relaxation equation given in equation (31) is 
then used to compute the new vorticity value. The stream function is then found by 
computing 


Sn - Sn-1 + 



m,n 


(78) 


from the wall to the edge of the boundary layer with the relaxed value of the stream func- 
tion computed from equation (33). 

The linear stability of the global iterative scheme for the vorticity-transport equa- 
tion is examined in appendix B for forward flow. For reversed flow the same scheme is 
used for the prescribed wall-shear calculations that is used for the prescribed displace- 
ment thickness calculations; this scheme is shown in appendix A to be stable for r S 1. 
For forward flow, the analyses are also quite similar, since for r S 1, the global 
iterative scheme is stable if the grid is refined to suppress the possible destabilizing 
influence of a lower order term. In practice, no instabilities which could be traced to a 
violation of the von Neumann condition were observed. 


Formulations 2 and 3 .- The numerical treatment of the vorticity-transport equation 
in these two formulations is the same as that used for formulation 1, with some minor 
changes. Small oscillations present in the solution of formulations 2 and 3 were elimi- 
nated by representing the term (or 8^/8x) by a combined first- and second- 

order difference expression which is given by 


di{/ 


= /3l 


^m,n ^ m- 


l,n ’^m- 


m,n 


2 A| 


^ (1 - 


- ij/. 


m- 


(79) 
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where /3 is a weighting factor in the rai^e of 0 = /3 = 1. If /3 = 1, equation (79) has 
second-order accuracy; if ^ = 0, the equation has first-order accuracy. The value of /3 
typically used was 0.8 or 0.9. 

In contrast to equations (44) and (52), the transformed shear in formulation 3 appears 
undifferentiated in equation (66). In the finite-difference representation of equation (66), 
this term was treated implicitly. For example, equation (75) which is for forward flow 
becomes 


Bjj + Q! — 2 + 3Cx + 2 Cy + Cg 


(80) 


where 


Cz = Atj^u 


f3m 


m 


m,nl 


(81) 


and the definitions of Cx and Cy are analogous to those expressions for and C-q 
given in equation (76). 

In both formulations 2 and 3 the terms in equations (55) and (70) are represented by 
central differences and the resulting equations are solved for u (or u) assuming co 
(or t) is known. The stream function is then deduced from equation (58) or (73) by using 
the trapezoidal rule. 

In order to complete the solution at a given ^ station in formulation 2, equa- 
tion (51) is solved to first order for Ug to give 


u^"*"^ = u^ 
e,m e,m 


+ r 


uj"^^ 1 

e,m-l 



- 

m ^3 


3cu 


m. 


(82) 


where j = o>o(4) which is the prescribed wall vorticity and r is the relaxation fac- 
tor. The solution changed only slightly when equation (82) was modified to give second- 
order accuracy. In formulation 3 the pressure-gradient parameter m is deduced in a 
similar manner from 


m^+1 = mi 

m m 


+ r 




- 

- 3t^) - m^ 

2 Ay 

\ m,2 

m,3 0/ m 


(83) 


In both equations (82) and (83) the value of the relaxation factor used was typically 0.6, 
the same as that value used in solving for o) (or t). 


26 



RESULTS AND DISCUSSION 


Computational Rate 

The computational rate on the CDC 6600 computer system for the results presented 
in this section is about 2300 grid points per second for the prescribed displacement 
thickness cases and about 3000 grid points per second for the prescribed wall-shear 
cases. In the global iterative procedure the number of iterations required for conver- 
gence is typically between 100 to 200; the prescribed initial conditions are the upstream 
boundary conditions repeated at each streamwise station. The computer time required 
for these calculations is usually 5 minutes or less, depending on the grid and the total 
number of iterations. In the prescribed displacement thickness procedure, the number of 
global iterations required for convergence was reduced considerably by using improved 
initial conditions obtained from the forward-marching procedure, which requires an aver- 
age of 20 to 30 column iterations at each streamwise station. 

Prescribed Displacement Thickness 

Comparison with solutions of the Navier - Stokes equations .- Comparisons with 
Briley's (ref. 29) numerical solutions of the Navier-Stokes equations for the flow in the 
vicinity of a laminar separation bubble are presented. Figure 2 shows a schematic dia- 
gram of the computational region used by Briley in solving the vorticity — stream-function 
formulation of the Navier-Stokes equations. The outer boundary was located about four 
upstream boundary -layer thicknesses from the wall. Along this boundary, a linearly 
retarded boundary velocity was prescribed, followed by a constant value. This 
boundary condition results in flow separation and reattachment as shown schematically 
in figure 2. Briley referred to this prescribed boundary velocity as \x^; however, in the 
present report u^ is the streamwise velocity component evaluated at the boundary -layer 
edge which is closer to the surface than the position of the outer boundary used by Briley. 
At the upstream boundary which was located at x = 0.05, Howarth's (ref. 30) solution of 
the boundary-layer equations was used for the boundary conditions. Briley discusses the 
fact that his solution of the Navier-Stokes equations does not account for the interaction of 
the viscous solution on the inviscid flow; likewise, the present boundary-layer approach 
has not been coupled to an inviscid calculation. Hence, both calculations can be regarded 
as the first stage of iteration with an inviscid flow. It is of interest to compare the re- 
sulting solutions in the viscous region. In the present calculation Howarth's solution is 
also used for the upstream boundary condition. The other boundary conditions are given 
in figure 1. 
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Figure 2.- Computational region used by Briley for Navier-Stokes equations. 


In figures 3 and 4, the present results are compared with those obtained by Briley 
(solution 4 in ref. 29); in general, the agreement between the two solutions is excellent. 
The reference Reynolds number which Briley used is Roo^l “ ^ Figure 3 shows 

the displacement thickness deduced by Briley which is prescribed for the present calcula- 
tions. A comparison of streamlines is given in figure 3 where = 0 is the divid- 

ing streamline which separates the inner recirculating flow from the outer main flow. 
Figure 4 shows good agreement in the skin friction obtained by the two approaches. The 
skin-friction coefficient is related to the shear by Cf ^R oo^l = The grid spacings 
used in the present solution in figures 3 and 4 are A| = 0.00625 and Ari = 0.05; only 
slight changes were observed in the solution when these values were doubled. On the fine 
grid, starting with crude initial conditions, the number of iterations required for conver- 
gence was 131; on the coarse grid, the number required was 106. 

Figure 5(a) shows the boundary -layer edge velocity Ug deduced in the present cal- 
culation and the outer boundary velocity uij prescribed by Briley. The agreement 
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Figure 4.- Comparison of skin friction. 
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(a) Comparison of boundary and edge velocity. 

Figure 5.- Velocity distributions and Meksyn pressure gradient parameter. 
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between the two distributions is relatively good. This good agreement suggests that the 
variation in the flow quantities between the boundary -layer edge and the outer boundary 
used by Briley is fairly small. However, according to Meksyn's (ref. 31) criterion, if Ug 
were exactly equal to U]-, then the present boundary -layer solution would be singular at 
the separation point. From his approximate asymptotic method for the boundary-layer 
equations as well as from the experimental data of Schubauer (ref. 32) and others, Meksyn 
concluded that a necessary condition for a regular solution of the boundary -layer equa- 
tions at separation is that the function 



j \ 

must have a local minimum prior to separation; hence, — - > 0 at Cf = 0. Figure 5(b) 

dx ^ 

shows that the present solution satisfies this requirement; if the edge velocity were the 
same as Briley's boundary velocity, then the resulting \ does not. All solutions in the 
present investigation were found to satisfy Meksyn's criterion at the separation point. 
Klineberg and Steger found the same result in their calculations of separated flow, al- 
though they based their conclusions on the pressure-gradient parameter 


m 


X ^^e 
^e dx 


(85) 


which is approximately equal to 

The calculation of the edge velocity serves as a check on the solution accuracy 
since, as noted in the prescribed wall-shear formulation, the edge velocity can be deduced 
either from the velocity profile or from the momentum equation evaluated at 77 = 0 re- 
sulting in 


d^e _ 

^ dx 5* a?7 


( 86 ) 


Knowing the wall- shear derivative and the value of u^ at the upstream boundary, equa- 
tion (86) ii^ integrated with respect to x (as given in eq. (57)) to obtain the Ue distribu- 
tion. In the present computation of the Briley case, approximately a 10-percent difference 
in the two values of Ug was found at each x station when a first-order scheme was 
used for the calculation in the forward flow region. This difference was reduced to about 
0.1 percent for the same grid when the second-order accurate Crank-Nicholson scheme 
was used. This result demonstrates the value of a second-order scheme in the forward 
flow region as it avoids the excessive number of grid points required to obtain an accurate 
solution with a first-order scheme. 
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Co mpa rison with other ^und ary-layer solutions .- In this section a comparison be- 
tween a solution obtained by Klineberg and Steger using a point iterative scheme for a 
prescribed wall shear and a solution using the present technique with 6* prescribed is 
presented. The transformed wall shear prescribed by Klineberg and Steger is related to 
the usual vorticity in equation (65) and is given by 

^o(x) = - 2)(X - 6) = Tj(x) (87) 

except for 2 = x = 6 where 

1 + 0 !(x - 2)(x - 6) (88) 

where a is a given parameter. For o? = 0.1, figure 6 shows the displacement thickness 
which they deduced for a typical separated flow calculation. The Blasius flat-plate dis- 
placement thickness is also shown in figure 6 for a basis of reference. If Klineberg and 
Steger's 6*(x) is used as a prescribed condition and if their solution at x = 1 for the 
upstream boundary conditions is used the skin friction deduced in the present calculation 
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Figure 6.- Displacement thickness distributions. 
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Figure 7.- Comparison of skin friction and pressure-gradient parameter. 

is shown to be the same as that prescribed by KLineberg and Steger (fig. 7). Excellent 
agreement is also shown in figure 7 in comparing the pressure-gradient parameter 
(x/ue)due/dx obtained in the present solution with that found by Klineberg and Steger. 

The corresponding values of Ug agree to a few tenths of 1 percent. Profiles of the u 
component of velocity are plotted in figure 8 at several streamwise stations. No differ- 
ence could be distinguished between these profiles and those obtained by Klineberg and 
Steger. Thus, at least for this case of separated flow, two inverse boundary-layer solu- 
tion techniques which are significantly different produce essentially the same solution. 

Using the same convergence criterion, this case of separated flow required 216 iter- 
ations for convergence in the present calculations as compared to approximately 850 re- 
quired for convergence by Klineberg and Steger. This difference may be caused by the 
small underrelaxation factor (0.05) used by Klineberg and Steger on the pressure-gradient 
parameter in the early stages of iteration. No underrelaxation was needed for this case in 
the present calculations. The unknown edge velocity and the pressure gradient are elimi- 
nated from the present formulation; such elimination should increase the convergence rate 
caused by the sensitivity of the boundary-layer solution to the pressure gradient. In addi- 
tion, the point iterative scheme may be slower since there is a delay in satisfying either 
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Figure 8.- Velocity profiles at several streamwise stations 
from present solution (6* prescribed). 

the wall or the edge boundary condition, depending on the direction in which the solution is 
updated. In the column iterative scheme these boundary conditions are imposed simulta- 
neously in updating the solution at each streamwise station. The above difference in the 
total number of iterations is partly compensated by the fact that the point iterative scheme 
requires less computer time per grid point than the column iterative procedure. 

Adchtional calculations.- As a further test of the prescribed 5* technique, addi- 
tional calculations were performed in which the approach to separation was much more 
rapid and the reverse flow velocities were larger than in the Briley and Klineberg and 
Steger calculations. The prescribed 5* used in these additional calculations for ^ 

= is given by the cubic equation 


^*(1) - ” ^o) + ~ ^o) ' ^o)^ 


(89) 


The constants a^^ to a^ are determined so that, at ^ = ^q, the prescribed 6* distri- 
bution matches the value and slope of the Blasius flat- plate distribution and, at ^ 
reaches a maximum value. For ^ given by a second cubic equation 


6*(^) = 5i + R2(^ - ^i) + a3(| - + ^4(4 - ^4)^ (90) 

where a^ to a.4 are determined by matching the first cubic at | = At ^ 
becomes a minimum value 62*, which is generally chosen about the same as the Blasius 
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Figure 9.- Prescribed displacement thickness distributions. 


flat-plate value at that particular ^ location. In figure 9, two 6* distributions given 
by equations (89) and (90) are shown; in these equations, = 1.065, = 1.35, and 

^2 = 1.884. The only difference between these two cases is that 1*^ case A 

and in case B. For these two cases (d6*/d^)j^ax “ 20 and (d5*/d^)jj^a^ 

= 36 as compared to the values of 12 and 4.3 in the Briley and Klineberg and Steger cal- 
culations, respectively. The displacement thickness adds an angle to the given body shape 
d 6^ X ^ 

by an amount . For example, if R„ t = 10^, then in case B at the point of 

(d6*/d^)j^^ the effective body slope is increased by 6.5°. For an aerodynamic shape 
such as an airfoil, a change in body slope of only a few degrees results in significant 
changes in the pressure distribution. 

Figure 10 shows the streamline pattern in the shear layer and recirculation region 
for case A. The thickness of the reversed flow region is amplified here since the normal 
coordinate has been multiplied by the square root of the Reynolds number. This separa- 
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Figure 10.- Computed streamline pattern in separation bubble for case A. 


tion bubble, which is typical of those calculated in the present investigation, is actually 
thin and is entirely confined to the interior of the viscous layer. 


Figures 11 and 12 show the resulting solutions for the skin-friction, edge velocity, 
and pressure-gradient parameter X for cases A and B, respectively. Only slight differ- 
ences are observed in these calculations for the two different grid spacings which were 
used; hence, the present technique appears to yield accurate solutions on a relatively 
coarse grid. Comparison of the skin-friction distributions in figures 11(a) and 12(a) 
shows that in case B, because of the greater 6*, the separated region is of greater extent 
than that in case A. The edge-velocity distributions for these two cases (figs. 11(b) and 
12(b)) illustrate the usual variations found in a flow involving a significant amount of sepa- 
ration. Near the points of separation and reattachment, the edge velocity or pressure 


coefficient 



varies rapidly with an interim region of relatively constant 


pressure which is usually referred to as the pressure plateau region. The plateau region 
is better defined in case B since the separated region is larger. 


In figures 11(c) and 12(c) it is seen that the pressure-gradient parameter X given 
in equation (84) has a positive slope at separation and, hence, satisfies Meksyn's criterion 
for a regular solution at the separation point. In figures 11(c) and 12(c), as well as in 
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figures 5(b) and 7 for the Briley and Klineberg and Steger calculations, respectively, a 
second local minimum occurs in X (or m) just prior to reattachment. Similarly, in a 
recent investigation, Horton (ref. 13) presented a calculation for a prescribed wall shear 
in which Meksyn's criterion is satisfied at both separation and reattachment. These cal- 
culations provide further support for Meksyn's observation that the shape of the X dis- 
tribution is similar for a wide range of separated flows. 
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(b) Edge velocity. 

Figure 11.- Effect of grid size on calculations for case A. 










X 


(c) Meksyn pressure-gradient parameter (S = Separation, R = Reattachment). 

Figure 12.- Concluded. 

Approximate solutions .- All of the calculations discussed herein were recomputed 
using the approximate forward- marching procedure. In all cases the calculations were 
stable and, as expected, the agreement between these approximate solutions and the more 
exact calculations obtained with global iteration becomes poorer as the magnitude of the 
reverse flow velocity increases. The solutions were almost identical in the Briley and 
Klineberg and Steger cases. The magnitude of the skin-friction distributions for cases A 
and B obtained by the forward- marching technique and by the global iterative procedure 
are shown in figure 13. The agreement is quite good in case A where the maximum re- 
verse flow velocity is -O.OSUoo; in case B where the velocity is about -O.lOUoo, the accu- 
racy of the reversed flow approximation is slightly less. For case B the approximate 
technique predicts about the same reattachment point as the more accurate global scheme 
(fig. 13). Therefore, at least for this case, the solution beyond reattachment does not 
appear to be affected by the local errors made in the reversed flow region. 
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4- Global iteration 


Forward-marching procedure 
u<0,set = 0 



Figure 13.- Comparison of skin friction for global and forward- marching 

procedure for cases A and B. 


Klineberg and Steger also converted their global iterative procedure in which the 
wall shear is prescribed to a forward-marching technique. They included the Reyhner 
and Fliigge-Lotz approximation and used a backward difference on 8u/9x in the continu- 
ity equation. However, their calculations were unstable except when a coarse grid was 
used with a first-order scheme in the stream direction. The instability occurred even 
though the magnitude of the reverse flow velocity was less than 0.02Uoo. In contrast, in 
the present calculations, some of which are for more severe separation cases, no insta- 
bilities were encountered even with further grid refinements, when the streamwise con- 
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vection of vorticity in the reversed flow region was neglected. Since both the reversed 
flow approximations and the problem formulations are different in these two studies, it is 
difficult to assess precisely why the present results are stable and those of Klineberg and 
Steger are not. 

The reduction in computer storage obtained by using the Reyhner and Fliigge-Lotz 
approximation is substantial since only two columns of data need to be stored. In con- 
trast, an entire plane of data is required in the global procedure. Significant reductions 
in computer time are also achieved. For example, in the present calculation of Briley's 
case, the total computer time required in the approximate calculation for 5700 grid points 
was about 1 minute on the CDC 6600 computer system; 5.5 minutes (131 iterations) were 
required in the global iterative procedure. Both of these computer times are significantly 
less than the 45 minutes used by Briley for the Navier-Stokes equations for 1050 grid 
points on the UNIVAC 1108. The computational speed of the CDC 6600 is approximately 
four times that of the UNIVAC 1108. 

It was found that the convergence of the global procedure could be accelerated con- 
siderably by using the approximate solution as initial conditions. Table I shows the num- 
ber of global iterations required for convergence when crude initial conditions (upstream 
boundary condition repeated at each streamwise station) are used in comparison with the 
use of the approximate solution. The average number of column iterations required at 
each X station in the forward- marching procedure are also shown in table I. As ex- 
pected, table I shows that this convergence acceleration becomes less as the magnitude of 
the reverse flow velocity increases and the reversed flow approximation becomes poorer. 
Use of the crude initial conditions resulted in an unstable calculation in case B; whereas, 
a converged solution was found by using improved initial conditions. 

TABLE I.- EFFECT OF INITIAL CONDITIONS ON TOTAL NUMBER OF 
GLOBAL ITERATIONS WITH 5 * PRESCRIBED 



Number of global iterations 

Average number of 

Maximum 

negative, 

u 

Case 

Crude initial 
conditions 

Initial conditions from, 
forward- marching 
procedure 

column iterations 
forward- marching 
procedure 

Briley 

131 

26 

10 

-0.01 

Klineberg and Steger 
(a = 0.1) 

216 

27 

13 

-.02 

A, coarse grid 

145 

29 

18 

-.05 

A, fine grid 

205 

114 

14 

-.05 

B, coarse grid 

Unstable 

130 

41 

-.10 

B, fine grid 

Unstable 

166 

j 

28 

-.10 
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Prescribed Wall Shear 


Comparison with solutions of the Navier-Stokes equations .- The wall shear computed 
by Briley (fig. 4) is used as a prescribed boundary condition for formulation 1. The re- 
sulting displacement thickness is shown in figure 14 to be in excellent agreement with that 
obtained by Briley. Similar agreement of the velocity profiles at various streamwise sta- 
tions is also found. This prescribed wall-shear calculation for the Briley case required 
219 iterations for convergence as compared to 106 and 131 iterations which were required 
for the coarse and fine grid calculations, respectively, with 6* prescribed. This result 
is somewhat surprising as it was originally thought that the solution would converge more 
rapidly with the wall shear prescribed than with 6* prescribed, since in the latter for- 
mulation the wall vorticity is unknown and must be deduced along with the solution. 


2.4 I — 


2.0 


1.6 


1.2 


6 * 


.8 


.4 



0 .1 .2 .3 .4 .5 

x/L 

Figure 14.- Comparison of computed displacement thickness distributions. 


Discontinuous solutions; formulation 1 .- Despite the excellent agreement obtained 
in the comparison with Briley's calculation, in several additional calculations with formu- 
lation 1, solutions were obtained which contain a streamwise discontinuity. For example, 
the previously discussed Klineberg and Steger case (for a = 0. 1) was calculated with the 
skin-friction distribution shown in figure 7 as a prescribed condition. In figure 15 the re- 
sulting solution for the edge velocity Ug obtained from the velocity profile is compared 
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Klineberg and Steger solution 

( 1^0 = 1/Ug^x/Ug j O)^ prescribed) 

Present solution 
(6* prescribed) 


1.2 


+ Present solution 

(co^ prescribed, formulation 1) 
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x/L 
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Figure 15. - Discontinuous solution for edge velocity, o? = 0.1. 
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with that found by Klineberg and Steger; their solution is the same as that deduced in the 
present formulation with 6* prescribed. A discontinuity in Ug is observed at x = 2.4, 
a point which is a short distance downstream of the prescribed separation point. The edge 
velocity Uq is negative between the discontinuity and the prescribed reattachment point 
X == 6.0. In this region the u component of velocity varies monotonically from the wall 
to the asymptotic edge value shown in figure 15; hence, in this region the entire boundary- 
layer flow is in the negative x direction. At the x location of the discontinuity, the v 
component of velocity has a sudden increase followed by a rapid decrease. This large out- 
flow from the boundary layer in the vicinity of the discontinuity is expected, since most of 
the flow is forward just before the discontinuity and all the flow is reversed just after the 
discontinuity. It should be noted that this calculation required 609 iterations to converge. 
The same results were obtained when the calculation was repeated with Ax increased 
from 0. 1 to 0.2. 

This calculation converged. However, it is not a valid solution to the boundary-layer 
equations, since the profile value of Ug (obtained from eq. (50) with 77 ^ °o) does not 
agree with that obtained from evaluating the x momentum equation at the wall (given in 
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eq. (51)). In formulation 2 these two values of Ue are forced to be equal. Using this 
formulation, this case was recomputed with the resulting Ug distribution smooth and in 
excellent agreement with that obtained by Klineberg and Steger. Good agreement is ob- 
tained in the resulting 5* distributions which are shown in figure 16. The present cal- 
culation required about 500 iterations to converge. As noted in the Briley case already 
discussed, this number of iterations is considerably more than the 216 iterations required 
in the prescribed 6* calculation. Thus, it is deduced that, in certain calculations, spec- 
ification of zero vorticity at the boundary -layer edge is not an adequate boundary condition 
to obtain a valid solution. Rather, it appears necessary to require the unknown velocity at 
the boundary -layer edge to be the same as that deduced from the pressure gradient in 
equation (51). In contrast, no discontinuities were encountered in any of the calculations 
with 6* prescribed since, in a similar manner, the flow in the boundary layer is con- 
trolled by setting the transformed stream function if/ equal to zero at the outer boundary. 



Figure 16.- Comparison of computed displacement thickness distributions, a = 0.1. 
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Discontinuous solutions; formulation 3 . - Klineberg and Steger encountered an unex- 
plained discontinuity in the pressure- gradient parameter when they set a = 0 in the pre- 
scribed wall shear given in equations (87) and (88). This transformed shear distribution 
and that distribution for a = 0.1 used in the calculation discussed above are shown in 
figure 17. Both of these shear distributions were used as input conditions to the present 
solution technique using formulation 3. For a = 0.1, the calculation converged in about 
500 iterations and agrees well with the solution obtained by Klineberg and Steger. For 
Q! = 0, a discontinuity occurs in the m distribution which is located at approximately 
the same x position as that found by Klineberg and Steger (fig. 18). The present 
numerical technique and problem formulation differ significantly from that of Klineberg 
and Steger and yet both procedures deduce essentially the same discontinuous solution 
for o! = 0. Therefore, it is concluded from this calculation that an arbitrarily pre- 
scribed wall-shear distribution does not necessarily result in a physically meaningful 
solution. The cause of this discontinuity is not known. Klineberg and Steger concluded 
that the discontinuous solution shown in figure 18 may indicate a possible breakdown in the 
boundary -layer equations. The breakdown may be caused by the relatively large extent of 
separated flow and the magnitude of the normal component of velocity at the boundary- 
layer edge. However, this argument is discounted since several cases with smooth solu- 
tions involving much larger normal velocities were computed in the present investigations 
with 6* prescribed. 



Figure 17.- Prescribed transformed wall-shear distributions. 
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Klineberg and Steger solution 

(f Q prescribed) 

Present solution 

(Tq prescribed, formulation 3) 



CONCLUDING REMARKS 

This report demonstrates that regular solutions of the laminar boundary -layer equa- 
tions can be obtained at separation with either a prescribed displacement thickness or 
wall-shear distribution. Excellent agreement was obtained in a comparison of these solu- 
tion techniques with a calculation of the Navier-Stokes equations made by Briley for a flow 
having a relatively thin separation bubble terminated by reattachment on a solid surface. 
Based on this agreement, as well as on that obtained by other investigators for compress- 
ible flows, it is concluded that the boundary -layer approximations to the Navier-Stokes 
equations are adequate for separated laminar flows of limited extent. 

The present inverse boundary-layer solutions are obtained by several different 
global iteration schemes in which the finite-difference scheme is switched in the reversed 
flow region to account properly for the flow direction. These global schemes have uncon- 
ditional diagonal dominance and are shown by the von Neumann analysis to have linear 
stability. As a special case in this analysis, it is concluded that overrelaxation of the 
diffusion equation results in an unstable calculation, at least with the three finite - 
difference schemes considered. In the forward flow region, it was essential to use a 
second-order accurate technique; in the cases considered in this report, first-order 
accuracy was found adequate in the reversed flow region. 
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The prescribed displacement thickness method described herein is preferred to the 
three schemes presented for a prescribed wall shear. For the same case, the results 
obtained with the wall shear prescribed required approximately twice the number of itera- 
tions and twice the computer time as that required with the displacement thickness pre- 
scribed. Furthermore, in the wall-shear prescribed formulation, it was necessary to 
constrain the unknown edge velocity so as to match that velocity deduced from the unknown 
pressure gradient. Even with this constraint, a discontinuous solution was obtained with 
the present formulation. The discontinuous solution was found to be similar to that ob- 
tained by Klineberg and Steger for the same case. Since the formulations and numerical 
techniques differ significantly for these two wall-shear prescribed procedures, it is con- 
cluded that an arbitrarily prescribed wall shear distribution does not necessarily result 
in a physically meaningful solution in separated flow. Further research is needed to 
determine the cause of this streamwise discontinuous behavior which occurs just down- 
stream of the separation point. 

An approximate forward- marching procedure is presented. In this procedure, a 
separated boundary-layer flow is computed in one sweep from the upstream boundary, that 
is, in the same manner as an attached flow. In the reversed flow region, the streamwise 
convection of vorticity is set equal to zero to prevent the instability which occurs in march- 
ing against the flow. This approximation is similar to that, first introduced by Reyhner and 
Fliigge-Lotz for reversed flow, of neglecting the streamwise convection of momentum. In 
this approximate scheme, the usual column iterative procedure was modified to give un- 
conditional diagonal dominance. The resulting column iterative scheme is analyzed for 
linear stability by the von Neumann analysis. This procedure eliminates the need for 
artificial convection terms which Reyhner and Flugge-Lotz and Werle, Polak, and Bertke 
added to obtain convergence in their column iterative procedure at each streamwise 
station. 

Comparisons between results obtained with the global iterative procedure and the 
forward-marching procedure show that the latter is accurate provided the magnitude of 
the reverse flow velocity is less than about O.lOUoo. The local errors made in the re- 
versed flow region with this approximate scheme are not observed to affect the flow down- 
stream of reattachment, at least for the cases considered. As expected, the number of 
iterations required for convergence in the global procedure is reduced substantially by 
using the approximate solution as initial conditions. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
September 5, 1975 
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APPENDIX A 


STABILITY ANALYSIS OF PRESCRIBED DISPLACEMENT THICKNESS METHOD 

Attached Flow 

The von Neumann stability analysis is applied to the global iterative scheme in 
equation (37) for the prescribed displacement thickness method for both attached and 
separated flows. Attached flow is considered first. Substituting for An, Bn, Cn, and 
Dn from equation (18) into equation (37) gives 


- 1 H- * i(c, - 


= oi(jO'^ + r 
m,n 


+ (1 - r) 




(Al) 


where the coefficients C^, C^j, and a, which are held constant in the von Neumann 
analysis, are evaluated with the most recent solution. A Fourier component of arbitrary 
wave numbers, and k-q, is written as 


U) 


j 

m,n 



(A2) 


where is the amplitude coefficient at the jth iteration. Substitution of equation (A2) 

into equation (Al) gives the amplification factor |g| = )N|y|D| where 


N = Q! + (1 - r)^l - cos 4> + 2C^ + i(l - r)C^ sin 4> 


D = O' + (1 + r cos 0)(1 - cos 0) + 2C^(1 - r cos 0) + rCq sin 0 sin 0> (A 3 ) 


+ 1 


r sin 0^2C^ - 1 + cos + Cq sin 0(1 + r cos 0) 


where 0 = k^ and 0 = k^ Ari. Note that the range of 0 and 0 values is 
0 = (0,0) = 7T and corresponds to the range of solution wavelengths and of 
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CO g 2 2 A4 




(A4) 


At7| 

since and can be expressed as 

k. 

k 

k?7 - — 

^71 


(A5) 


The requirement for stability is that | g | = 1. It is convenient to consider the stability 
requirement in the form Igp = 1 which becomes 


d2 - n 2 g 0 

Expressing equation (A3) in the form of equation (A6) results in 

r^l + cos 6)(1 - cos 4>)(a + 1 - cos <p) + sin 0\|l + cos 6 + 2C^\jl - cos 0 
+ 2C| Q!(l - cos 0) + 2(1 - r)(l - cos 0)J + sin 0 sin S 0 


(A6) 


(A7) 


Equation (A7) must be satisfied to prevent instability. 

Equation (A7) is first considered in the limit of and A77 ~ 0 with 

g+2 ^ 2 

held constant so as to preserve the highest ^ and 77 derivatives. Since 
a ~ Cfj ~ Atj then with a and set equal to zero, equation (A7) reduces 

(1 + cos 6)(1 - cos (p)^ + 4C^^(1 - cos 6) + 4C^(1 - r)(l - cos 4>) = 0 (A8) 

Equation (A8) must be satisfied for stability in the relaxation solution of the diffusion 
equation 


ug*2 8(4) _ 8 (4) 

" 87,2 


(A9) 
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and hence r = 1 for stability. Thus, overrelaxation cannot be used in obtainii^ a relaxa- 
tion solution of the diffusion equation with the Crank-Nicholson scheme. 

Equation (A7) is now examined to see if the lower order terms may result in the 
violation of this inequality. With r S 1 , and since and a are positive, then all of 
the terms in equation (A7) are positive with the exception of the last which is negative 
when Cjj < 0. The largest negative value of occurs as the boundary -layer edge is 
approached. In such a case, 

C77 - 1)^ (A12) 

where C-q < 0 if d 6 */d^ > 0. Therefore, as the positive slope of the prescribed dis- 
placement thickness increases, the destabilizing influence of the lower order term is in- 
creased. Obviously the effect of this term can be reduced by decreasing At]. In practice, 
the calculations were performed without regard to this observation; hence, the only sta- 
bility restriction is r = 1 which is deduced from the consideration of the higher order 
terms. This condition is consistent with the observation of Richtmyer and Morton that 
the stability of many finite -difference solutions of the diffusion equation is essentially 
unaffected by the presence of lower order terms. 

The only evidence of instability caused by the lower order terms occurred in the 
calculation in which the prescribed displacement thickness is case B shown in figure 9. 

This case is the most severe separation calculation considered in this report. In this 
case, very small amplitude oscillations occur in the 77 distribution of vorticity near the 
boundary-layer edge. In this region the negative values of Cq are sufficiently large to 
cause I g I to exceed slightly unity for certain values of <p and 6. The stability of the 
calculation is presumed to be enhanced by the imposition of the outer boundary conditions; 
of course, these conditions are not accounted for in the von Neumann analysis. It was 
observed that the magnitude of these oscillations was reduced by decreasing A 77 as 
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anticipated from equation (A7). Furthermore, the oscillations were entirely eliminated 
by using a first-order scheme in which a von Neumann analysis shows that the lower 
order terms do not have any destabilizing effect. 

The same calculation discussed in the preceding paragraph was used to test several 
values of the relaxation factor. With reasonably accurate initial conditions provided by an 
approximate forward-marching procedure, calculations were made for the different values 
of the relaxation factor of 0.4, 7, 1.0, 1.1, and 1.5. Only in the case of r = 0.4 was con- 
vergence observed. However, the behavior of the maximum residual (increment in vor- 
ticity between two successive iterations) is distinctly different depending on whether r = 1 
or r > 1. In both cases the maximum vorticity residual occurs at the surface. However, 
for 0.4 < r = 1.0 this residual, although small, has an oscillatory behavior and conver- 
gence is never obtained. In contrast, for r > 1.0 the solution diverges rapidly in a man- 
ner which suggests that the von Neumann stability criterion has been violated. It is con- 
cluded from these calculations that whereas the von Neumann analysis requires r = 1 as 
a necessary condition for stability, the actual value of r may need to be less than unity 
because of the calculation of the surface vorticity. 


Alternate Scheme for Attached Flow 


For r = 1, equation (Al) can be interpreted as the finite-difference representation 
of the unsteady vorticity-transport equation 


8u) 8u) , 8u) _ 8^0; 

8t ^ 8^ ^ Bn~'^ 


(A13) 


where the time derivative is 


Boj 

8t 




j + 1 - J 


At 


(A14) 


with At = Arj^/a and where a and b represent the coefficients in equation (12). In 
an attempt to eliminate the destabilizing effect of the lower order terms found in the pre- 
vious scheme the time derivative, which is introduced to establish diagonal dominance, is 
given by 


B(jj 

~Bt 


' CO 


- J 


m,n 


At 


m,n ^ 


CO 


- 


i+1 

m-l,n '~^m-l,n 
At 


(A15) 


A comparison of the computational molecules of these two schemes in figure 19 shows the 
symmetry of the alternate scheme. This alternate global scheme can be written as fol- 
lows by using equation (A15) with a = Atj^/2 At and after introducing a relaxation factor 
as was done previously: 
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u > 0 

Forward flow 

(alternate) 

Figure 19.- Time-like interpretation of iterative finite-difference schemes. 




au)^ + raco^ 


j -r ± LZ \ju ^ 

m,n m-l,n 


+ r 


-{Crj + ^ - 1 ~ 


+ (1 - r) 


!/■<-, . ^ + (o.n. . ^ 


(A16) 


i S lXn-1 - 2C^ . lj<,n 


The amplification factor for this scheme is | g | = | N |^]D| where 

N = a(l + r cos e) + (1 - r)(2C^ + 1 - cos + ij^-ar sin 6 + (1 - T)Cf] sin 

D = a(l + r cos 6) + (1 + r cos 0)(1 - cos cp) + 2C^(1 - r cos 6) + rCy sin <p sin (A17) 
+ ij"r sin e(2C^ - 1 + cos 0 - + C^j sin <p{l + r cos 6) 


The requirement for stability jgp = 1 becomes 
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= r|(l + cos 0)(1 - cos 0)[(1 + r)o! + 1 - cos <^>] + (c-q sin ^\Jl + cos 6 
+ 2C|^1 - cos 6J + 2(1 - r)C^j^(l - cos 9)a + 2(1 - cos 0)J 


+ (1 - r)C-qa sin 6 sin ^ = 0 


(A18) 


In the limit of zero mesh size with At; 2/A4 held constant, this scheme reduces to the 
same as that previously considered and, hence, r = 1 is required for stability. With 
r = 1 all of the terms in equation (A18) are positive except for the last which is negative 
if r 1 and Cjj < 0. Thus, in this global iterative scheme the possible destabilizing 
influence of the lower order term is eliminated provided no underrelaxation is used. Un- 
fortunately, the advantage of this scheme could not be used in the present study since a 
smaller value of the relaxation factor was found to be necessary in the calculations, pri- 
marily because of the computation of the surface vorticity. 


Separated Flow 

The stability analysis of the global iterative scheme for separated flow proceeds in 
an analogous manner to that for attached flow. Substituting for A^, Bn, Cn, and Dn 
from equation (21) into equation (37) gives 



+ ^q; + 2 - 


2C 




O) 


j+1 

m,n+l 



- 


= + (1 - r) 

m,n ^ 

■(Cj? + l)^m,n-l + ■ ^^)‘^in,n + (^7 ' ‘ 

'■)^m,n+l 


_ 



- 2rC|0) 


J 

m+l,n 


(A19) 


Since < 0 for separated flow, it is convenient to set The substitution 

of equation (A2) into equation (A19) gives the amplification factor 







1 - r(l - cos 6) 

|C| 1 + (1 - r)(l - cos <^) + i 

(1 - r)Cjj sin p + rjc^j sin Q 


I- 


+ 1 - cos (p + iC^j sin p 



(A20) 
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Expressing equation (A20) in the form of equation (A6) results in 


r|(l - r)(l - cos 0)jc^p - (1 - r) sin <p sin + 




X (1 - cos (p)^ + ^ 


(1 - cos (p) + ■ cos 6) + jC||(l - cos <p)^ - r - (1 - r cos 0)J| = 0 

(A21) 


which must be satisfied to prevent instability. First the limit of zero mesh size is con- 
sidered, which again implies the solution of the diffusion equation by relaxation. Setting 
At/ and A| equal to zero with C| held constant results in 

(1 - r)(l - cos 6)|C^ j^+ 2 ~ ■ cos (1 - cos (p)\2 - r - (1 - r cos 0)j = 0 


(A22) 

Equation (A22) is satisfied if r S 1. With this restriction on the relaxation factor, all of 
the terms in equation (A21) are positive with the exception of the second term which is 
negative when Cjj > 0. However, it is shown here that the first three terms of equa- 
tion (21) are always positive; hence, if r S 1, then equation (A21) is always satisfied. 

The first three terms of equation (A21) can be rewritten as 

-,1/2 


- 


- 

2 


- 

|/(1 - r)(l - cos e)jc^ 


I2 - r . _ 

1 2 

1 

2 

(1 - r)p ^ ^^(1 - cos 0) 


- (1 - r) sin 6 sin cp 


Ct C 


J 


{A23) 


With r S 1 and > 0 (worst case), the quantity in braces in equation (A23) is positive 
and the von Neumann condition is satisfied. Thus, the first-order global iteration scheme 
used for separated flow is stable provided that overrelaxation, that is, when r > 1, is not 
used. 
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STABILITY ANALYSIS OF PRESCRIBED WALL- SHEAR METHOD 

The von Neumann stability analysis is applied in this appendix to the global iterative 
scheme used for the vorticity transport equation in the prescribed wall-shear method. 

Only attached flow is considered since the reversed flow scheme is the same as that used 
in the prescribed displacement thickness method. Substituting for An, Bn, Cn, and Dn 
from equation (75) into equation (37) gives 

-(C, * + (2 + 3C| + + (C„ - 

* - <-2,n) (Bl) 

where the coefficients C^, Crj, and a are given in equations (75) and (76). The sub- 
stitution of equation (A2) into equation (Bl) gives the amplification factor 


= n + (1 - r) 



a + 2(1 - r)(l - cos 0) + 3(1 - r)C^ + i|^2(l - r)Cn sin 0 


a + 2(1 - cos 0) + Ce 3 - r(4 cos 6 - cos 26) 

L j 

+ i 

1 

2Crj sin 0 + rC|(4 sin 6 - sin 26) 


(B2) 

In the case of no underrelaxation, that is, when r = 1, equation (B2) becomes 


a 


a + 2(1 - cos <p) + 2C^(1 - cos 6) + i 2 Ctj sin 0 + C^(4 sin 6 - sin 26) 


(B3) 


which satisfies the von Neumann condition [g | ^ 1 since = 0 and a = 0. 

For a more general value of r it is convenient to rewrite the von Neumann con- 
dition |g| = |N|/|D| as 
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g 0 

Expressing equation (B2) in the form of equation (B4) gives 


(B4) 


3-3 cos 6(2 - cos 6) + 2r(l - cos + 2 sin 6 sin 0(2 - cos 6)C^C fj 


+ (2 - r) sin^ + (1 - cos 0)j^3(2 - r) - 


4 cos 6 + cos 29 


Q!|^(l - cos 6)^C^ + (1 - cos 0) 


+ (2 - r)(l - cos 0)^ Mo 


(B5) 


In the limit of A| and A77 equal to zero with held constant, equation (B5) 
reduces to 


r{ 3 - 3 cos 0(2 - cos 0) + 2r(l - cos 0) 


C| + (1 - cos 0 ) 


3(2 - r) - 4 cos 0 + cos 20 


+ (2 - r)(l - cos 0)^ > = 0 


(B6) 


The most restrictive condition on r is deduced from equation (B6) in the limiting case 
of C^~0~0— 0 which results in r § 1. This result is the same as that found previ- 
ously in this report for several other relaxation schemes for the diffusion equation. 

Examination of equation (B5) with r < 1 shows that all of the terms are positive 
with the exception of the second term which is negative when < 0. Numerical evalua- 
tion of equation (B5) for C-^ < 0 and various values of the remaining parameters showed 
that this term can be sufficiently large to violate this inequality for r < 1. However, this 
destabilizing term is a lower order term and, hence, its magnitude can be reduced by 
decreasing A77. In practice, no instabilities caused by this lower order term were 
encountered in the calculations. 
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